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Grid-Free  Electrostatic  Particle  Simulations 

Introduction  This  report  summarizes  work  done  between  September  2007  and  June  2008  re¬ 
lated  to  the  development  of  a  set  of  grid-free  simulation  tools  for  problems  in  plasma  dynamics. 
As  noted  in  the  abstract  on  page  one,  the  research  objectives  are:  (1)  to  develop  a  grid-free  field 
solver,  fluid  model,  and  kinetic  model,  (2)  to  evaluate  these  tools  in  comparison  with  traditional 
mesh-based  methods,  and  (3)  to  demonstrate  the  capability  of  the  grid-free  approach  in  an  appli¬ 
cation  of  USAF  interest. 

The  past  year’s  work  has  focused  on:  (1)  Grid-free  3D  Direct  Simulation  Monte  Carlo  and  opti¬ 
mal  parallel  implementation,  (2)  application  of  grid-free  Direct  Simulation  Monte  Carlo  methods 
to  the  continuous  wave  formation  of  a  Bose  Einstein  condensate  and  Hypersonic  Flows  and,  (3) 
Validation  of  grid-free  agents  MONICO  and  G2.  (4)  Boundary  integral  methods  and  Regularized 
Particles,  (5)  Boundary  Integral  Corrected  PIC-An  Advanced  Subcell  Method  (6)  Extended  the 
treecode  to  kernels  of  the  from  re~cr ,  ideal  for  spacecraft  plume  simulations.  The  work  accom¬ 
plished  has  provided  several  key  results,  as  described  below. 

Basic  plasma  science  plays  an  increasingly  significant  role  in  applications  of  importance  to  the 
United  States  Air  Force  (USAF).  These  include  the  study  of  satellites  and  how  they  interact  with 
an  ion  engine,  the  use  of  electronic  tethers  to  protect  satellites  and  how  solar  eruptions  effect  satel¬ 
lites.  Additional  examples  include  the  basic  physics  of  magnetrons  for  use  in  directed  energy  (DE) 
systems,  the  prediction  and  prevention  of  arc  formation  in  DE  systems,  understanding  intense  laser- 
plasma  interaction  as  it  relates  to  space  propulsion,  “non-destructive”  materials  testing  for  turbine 
blades  in  aircraft,  and  the  study  of  Penning  traps  as  a  potential  way  of  trapping  anti-protons  and 
positrons  to  create  and  store  anti-hydrogen  as  a  fuel,  etc. 

The  above  plasma  applications  are  examples  of  plasmas  which  require  a  fully  kinetic  description  in 
at  least  part  of  the  domain.  The  most  common  approach  to  simulating  these  types  of  systems  is  to 
make  use  of  a  fully  Lagrangian  framework,  where  the  model  is  reduced  to  tracking  the  evolution  of 
test  particles  in  phase  space.  Of  the  many  varieties,  the  most  accepted  approach  is  Particle-In-Cell 
(PIC)  [1].  The  PI,  with  his  collaborators,  is  developing  an  alternative  approach,  the  Boundary  In¬ 
tegral  Treecode  (BIT)  [2]  and  Gridless  DSMC.  BIT  is  based  on  fast  summation  algorithms,  which 
eliminates  the  need  for  a  fixed  volumetric  mesh.  Gridless  DSMC  uses  the  clustering  strategy  of 
BIT  to  determine  collision  partners  and  combines  this  with  mesh  free  interpolation  to  track  flow 
averages  [3,4]. 

In  a  range  of  numerical  experiments,  the  all  scale  resolution  of  BIT  has  proven  to  provide  a  sub¬ 
stantial  improvement  over  traditional  PIC  [2].  Further,  the  Gridless  DSMC  has  demonstrated  the 
ability  to  more  accurately  capture  flow  fields  for  nonuniform  flow  fields,  such  as  in  hyper  sonic 
flows  [3,4].  The  Parallel  version  of  the  Gridless  DSMC  has  been  tested  for  scalability  and  is  cur¬ 
rently  scaling  at  about  90  percent  up-to  128  nodes  of  a  distributed  memory  machine,  and  it  has 
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not  been  test  beyond  this.  Further,  the  gridless  DSMC  appears  to  avoid  issues  of  cell  interface 
temperature  jumps  due  to  the  stochastic  nature  of  the  tree  clustering  algorithm.  The  randomness 
can  be  attributed  to  the  fact  that  the  clustering  is  combined  with  Monte  Carlo  collisions  which  adds 
randomness  to  particle  properties. 

The  BIT  method  has  been  shown  to  have  significantly  less  numerical  heating  than  explicit  PIC 
codes  [5].  The  objective  in  this  currant  work  is  to  develop  a  grid-free  electromagnetic  formulation 
of  BIT.  BIT  is  a  numerical  Green’s  function  method  based  on  discretization  of  the  integral  for¬ 
mulation  of  the  problem.  In  the  electromagnetic  case,  the  Green’s  function  is  time  dependent  and 
therefore  discretization  in  time  requires  that  a  time  history  be  maintained.  The  issue  of  a  time  his¬ 
tory  has  already  been  addressed  for  time  dependent  Green’s  functions  with  stationary  point  sources 
[7,8].  However,  for  moving  point  sources,  the  time  history  is  still  a  computational  issue,  since  ev¬ 
ery  step  introduces  N  new  sources  which  need  to  be  maintained  for  all  future  time.  One  approach 
to  handle  this  issue  is  the  to  use  a  coarse  fixed  underlying  mesh  for  tracking  the  effect  of  distant 
moving  charges. 

As  a  first  step  towards  the  electromagnetic  case,  we  are  exploring  the  use  of  a  coarse  fixed  mesh  for 
domain  decomposition  in  electrostatic  problems,  where  BIT  is  used  as  a  sub-cell  method  within 
each  PIC  cell  [6].  The  methods  presented  here  differ  from  Particle-Particle-Particle-Mesh  (P3M) 
in  that  local  boundary  integrals  are  used  within  each  cell  to  provide  an  accurate  description  of  the 
local  fields  within  a  mesh  cell  [9].  We  demonstrate  that,  in  ID,  regularized  BIT  corrected  PIC 
substantially  reduces  numerical  heating,  even  when  Ax  A>  A./>  We  expect  that  Ax  A>  Xd  should 
increase  the  efficiency  of  large  scale  PIC  calculations,  enhancing  the  capability  of  legacy  PIC 
codes  by  adding  both  fidelity  (less  numerical  heating)  and  efficiency  (allowing  more  flexibility  in 
meshing  at  a  given  level  of  fidelity).  For  a  wide  class  of  electromagnetic  problems,  where  PIC 
is  used,  this  electrostatic  work  is  beneficial,  since  for  stability,  the  mesh  spacing  is  sufficiently 
small  that  the  system  looks  electrostatic.  In  collaboration  with  AFRL/DEHE,  the  PI  is  working 
on  incorporating  this  into  ICEPIC.  This  work  is  intended  to  complement  our  planned  study  of  the 
use  of  fast  summation  methods  in  compressing  the  time  history  for  time  depended  moving  sources. 

New  Efforts  (Sep  2007  to  Present):  PI  Karsny  and  PI  Christlieb  have  been  working  on  the  extinc¬ 
tion  of  the  integral  methods  to  mesh  based  systems  for  the  ease  of  coupling  grid-free  solvers  with 
grid  based  solvers  and  we  are  working  on  three  papers.  PI  Karsny  has  put  substantial  effort  into 
the  development  of  recurrence  relations  for  kernels  of  the  from  re~cr  so  that  they  may  be  incor¬ 
porated  into  the  adaptive  treecode  framework,  which  he  is  in  the  process  of  writing  this  up.  This 
is  critical  for  the  modeling  of  neutral  as  well  as  dense  plasmas  simulations.  While  we  heave  not 
included  this  work  here,  it  will  stand  as  a  key  pice  of  a  white  paper  we  will  be  submitting  in  the 
near  future.  PI  Christlieb  has  focused  on  the  impact  of  regularization  on  solutions  in  BIT  and  BIT 
corrected  PIC  through  analysis,  development  of  high  order  C°°  ‘localized’  regularization  and  how 
it  improves  the  solution  in  BIT  and  BIT  corrected  PIC,  development  of  truly  localized  high  order 
(C'!)  regularization  in  BIT  (PI  and  Dr.  Ong),  the  devilment  and  analysis  of  high  order  correction 
inside  spectral  deferred  correction  (SDC)  (PI  and  Dr.  Qiu),  definition  of  a  new  metric  (accuracy 
criterion)  for  a  fair  comparison  between  single  step  methods  and  deferred  correction  methods  (PI, 
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Dr.  Ong  and  Dr.  Qiu),  the  application  12  order  SDC  to  BIT  simulations  of  one  and  two  stream  in¬ 
stability  (PI,  Dr.  Ong  and  Dr.  Qiu)  and  the  application  of  high  order  SDC  to  BIT  simulations  with 
point-insertion  (PI  and  Dr.  Ong).  Objective  of  the  work  with  Dr.  Ong  and  Dr.  Qiu  is  to  develop 
a  truly  high  order  particle  method.  Since  efforts  with  Dr.  Ong  and  Dr.  Qiu  are  now  just  coming 
into  there  own,  the  PI  will  details  these  in  the  2008  spring  review,  August  13-15,  for  computational 
math.  Additional  work  with  Dr.  Cartwright  has  focused  on  high  order  time  stepping  in  PIC  and 
the  possibility  of  truly  implementing  an  efficient  high  order  method  in  the  ICEPIC  framework. 
This  is  only  computationally  effect  because  of  the  information  that  is  already  stored  for  the  paral¬ 
lel  framework  and  the  thought  s  that  this  method  may  not  be  particularly  efficient  in  a  serial  PIC 
code.  Efforts  to  make  the  particle  push  not  have  a  self  force  is  the  primary  focus  of  this  effort.  The 
PI  has  been  collaborating  with  Dr.  Cartwright  on  this  effort  for  about  nine  months.  Finely,  the 
PI  and  Dr.  Olson  (a  former  student)  are  pursuing  the  application  of  Gridless  DSMC  to  hypersonic 
flows,  orbital  station  keeping  and  and  plume  modeling  of  trusters.  We  are  working  with  Dr.  Ong  on 
combining  the  3D  gridless  DSMC  with  3D  BIT  so  at  to  provide  accurate  mesh  free  simulations  of 
arc  formation.  On  another  note,  I  am  please  to  inform  AFOSR  that  my  former  student,  Dr.  Olson 
will  be  joining  AFRL  Kirtland  DEHE  as  a  full  time  research  scientist.  Dr.  Cartwright,  Dr.  Olson, 
Dr.  Ong  and  the  PI  are  starting  to  discuss  a  four  way  collaboration  on  arc  and  spark  formation  in 
wave  guides. 

Below  is  a  brief  outline  of  BIT  corrected  PIC.  Deities  of  the  work  can  be  found  in  an  attached 
paper  which  was  just  submitted  for  publication.  The  Gridless  DSMC  has  been  accepted  for 
publication.  Again,  a  brief  review  is  below.  Details  can  be  found  in  the  attached  preprint.  I 
have  also  include  a  biographical  sketch,  detailing  recent  accomplishments. 

Electrostatic  BIT  Given  a  collection  of  N  macro-particles  interacting  through  a  long  range 
coulomb  force,  the  question  becomes  what  is  the  most  efficient  way  to  evaluate  the  self  fields 
for  a  desired  accuracy.  For  a  collection  of  charged  particles  in  a  bounded  domain  the  electrostatic 
potential  is  given  by  the  solution  of, 

V2^(x)  =  -1^,  p(v)  =  £^-6(x-xy) ,  {x,x7'j  eQ/SQ  (1) 

to  7=1 

where  (|)  is  the  electrostatic  potential,  £0  is  the  permittivity  of  free  space,  p(x)  is  the  charge  density, 
{qj,xj}  is  the  charge/location  on  a  given  test  particle,  Q.  is  the  domain  and  r)Q  is  the  boundary  of 
Q.  The  problem  may  have  Dirclet,  Neumann,  mixed  or  periodic  boundary  conditions  along  d£l. 

In  general,  the  solution  to  equation  (1)  for  any  right  hand  side  may  be  expressed  as, 

4>(x)  =  -l  ^Y^G(x\y)dy  +  ^^(y)dnyG(x\y)+G(x\y)(dn^(y))ds^  (2) 

where  we  chosen  G  to  be  the  free  space  Green’s  function,  (])(y)  and  ( — 3„v(f)(y))  are  the  potential  and 
normal  electric  field  along  d£2.  The  problem  would  be  over  determined  if  both  (])(y)  and  (— r)„(|)(y)) 
were  specified  at  each  point  y  E  r)Q.  The  unknown  function  is  determined  by  imposing  consistence, 

lim  <|>(x)  =  <|>(y)  or  lim3%<|)(x)  =  3«v<|>(y)  (3) 

x^y  x^y  7 
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wherever  cf)(y)  or  (3„  d)(y))  is  specified  along  dQ. 

The  integral  equation  for  the  unknown  function,  y  or  a,  is  solved  using  numerical  quadrature  and 
the  resulting  matrix  is  inverted  using  a  modified  GMRES  [2].  The  resulting  system  looks  like  a 
collection  of  N  particles.  The  electric  field  at  particle  i  is  given  by, 

N  cr 

E(xj)  =  Y,  ^VxG(x|xj-)|x=xi  -  Vxfi.C.|x=Xi  (4) 

j=  hW  ° 

where  B.C.  stands  for  boundary  correction  term.  The  primary  issue  with  equation  (4)  is  that  the 
order  of  operations  count  is  0(N2).  In  most  simulations,  N  may  range  from  several  hundred  thou¬ 
sand  to  tens  of  millions  of  test  particles,  which  makes  the  calculation  prohibitively  expensive. 

To  overcome  the  0(N2)  issue,  BIT  makes  use  of  a  fast  summation  algorithm.  There  are  many  vari¬ 
ants  of  fast  summation  algorithms,  in  this  case  we  are  making  use  of  the  Treecode  algorithm.  The 
idea  is  to  approximate  the  long  range  interaction  of  each  distant  cluster  of  particles  as  a  moment 
expansions  about  the  center  of  each  cluster.  Since  plasma  problems  are  dynamic,  the  algorithm 
constructs  the  clusters  using  bisection  algorithm  recursively  in  space,  this  creates  a  hierarchal  tree 
sorting  of  the  particles.  Fields  are  computed  using  a  recursive  divide  and  conquer  approach,  which 
makes  use  of  the  tree  for  determining  the  appropriate  approximations  in  computing  the  long  range 
interactions  [2].  The  method  is  0(N\ogN),  a  substantial  speedup  over  direct  summation.  (Note 
that  PIC  is  0(M\ogM)  where  M  is  the  number  of  mesh  elements.) 

BIT  Corrected  PIC  Our  goal  is  to  provide  legacy  PIC  codes  with  a  robust  and  accurate  sub-cell 
method  for  field  calculations,  minimizing  numerical  heating  while  increasing  Av  A>  The  simple 
idea  is  to  use  regularized  BIT  as  a  sub-cell  method  within  each  PIC  cell.  The  boundary  conditions 
on  BIT  are  chosen  to  be  either  the  potential  or  field  computed  on  a  PIC  mesh.  The  domain  where 
BIT  is  applied  may  be  chosen  to  be  larger  than  a  single  cell.  Figure  1  a)  shows  how  the  force  falls 
off  between  two  particles  within  a  mesh  cell  (blue  dashed  line)  and  how  the  force  behaves  in  ID 
(solid  blue  line),  2D  (green  dashed  line)  and  in  3D  (red  solid  line).  Figure  1  b)  shows  a  cartoon 
of  a  PIC  mesh  where  the  red  cell  is  the  one  we  are  considering.  The  domain  where  BIT  is  applied 
could  be  as  small  as  the  red  cell  or  could  be  extended,  as  in  the  yellow  domain,  for  computing  the 
force  in  the  red  cell.  In  the  case  of  the  yellow  domain,  the  boundary  integral  uses  the  PIC  fields 
from  the  mesh  nodes  on  the  edge  of  the  domain  as  the  boundary  conditions  and  all  of  the  particles 
in  the  yellow  domain  are  used  in  computing  the  field  evaluations  for  the  particles  in  the  red  domain 
via  the  treecode. 


ID  Example  Consider  the  charge  neutral  test  problem  of  N  randomly  placed  negative  test  charges 
with  a  uniform  background  positive  charge  in  domain  Q.  =  ( a ,  b)  with  boundary  points  c)Q.  =  (a,  b}. 


d2§ 

dx2 


qo>iQ 

£o 


(b  —  a) 
N 


£8(x-x()j 


<K«)  =  <K*0  * 


d§(a) 

dx 


d§{b) 

dx 


Let  M  be  the  number  of  PIC  mesh  cells  and  let  the  mesh  spacing  be  Ay  =  (b  —  a)/(M  —  1).  Let 
the  mesh  points  be  give  by  y,-  =  z'Ay  4-  a,  i  e  (0, 1, . . .  —  1}.  Let  xj  be  the  location  of  the  test 
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Figure  1:  a)  is  a  plot  of  how  the  fields  fall  off  in  a  PIC  code,  as  a  function  of  distance  with  in 
the  cell,  as  compared  with  the  true  fields,  b)  is  a  cartoon  of  a  BIT  corrected  PIC  mesh  where 
the  red  cell  is  where  the  fields  are  computed  with  BIT  on  the  yellow  regain. 


particles  in  Q.  Let  (]),  be  the  mesh  based  solution  to  the  potential  at  mesh  cell  i  obtained  by  using 
an  FFT  (reference  [1])  and  let  C,  =  {xj\ [(xj  —  a)/ Ay \  =  /}.  The  sub-cell  electric  field  using  the 
regularized  Green’s  function  is  given  by, 


qim 
E0 

(b~a) 

2  N 


xk  Jj—  L  jeQ,  j^k 


Xk-Xj 


Zjsc,  vw 


-Xj)2  +  d2 


\/  ( xk-Xj)2+d 2 

\A 


:  +  d2 


(5) 


Figure  2:  a)  Force  Vs.  Location:  PIC  is  the  back  dashed  line,  BIT  is  the  green  line  with  dots,  BIT 
corrected  PIC  with  d  —  0  is  the  open  red  circles,  BIT  corrected  PIC  with  d  =  Av/10  is  the  blue 
dots  and  extend  cell  BIT  corrected  PIC  is  the  magenta  stars,  b)  FFT  of  force  for  three  results. 


Figure  2  a)  shows  the  solution  the  electric  field  computed  using  PIC  (the  black  dashed  line),  BIT 
(the  green  line  with  dots),  BIT  corrected  PIC  with  no  regularization  (the  open  red  circles),  BIT 
corrected  PIC  with  d  —  Ar/10  (the  blue  dots)  and  an  extended  cell  BIT  corrected  PIC  with  one 
cell  on  either  side  (magenta  stars).  The  simulation  uses  5  PIC  cells  and  N  =  500  randomly  placed 
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particles.  Observe  that  if  d  —  0,  BIT  corrected  PIC  is  in  excellent  agreement  with  BIT,  however,  if 
d  is  finite,  there  is  divergence  at  the  cell  interfaces.  The  divergence  at  the  cell  interfaces  is  because 
near  the  boundary,  the  regularization  pollutes  the  field  associated  with  charge  near  the  edge.  The 
mechanism  for  this  is  well  understood  by  the  PI,  the  extended  cell  represents  one  possible  solution. 
Alternatively,  instead  of  using  a  regularization,  one  could  propose  a  smooth  limiter  on  the  force, 
which  would  have  a  similar  effect.  Figure  2  b)  shows  the  FFT  of  the  force  in  ID  for  BIT  (identical 
to  direct  sum  in  ID)  BIT  corrected  PIC  with  extended  cells  and  regularization  and  Direct  Sum  with 
regularization.  Observe  that  the  spectrum  of  the  regularized  system  is  much  smoother. 


Figure  3  shows  the  results  for  several  numerical  heating  experiments  comparing  PIC,  BIT  and  BIT 
corrected  PIC  with  regularization.  The  system  is  charge  neutral  with  4000  electrons  uniformly 
sampled  in  space  and  the  velocity  is  samples  from  a  Gaussian  distribution  with  initial  temperature 
Te.  The  numerical  integrator  was  Leap  Frog  and  the  runs  are  100  plasma  oscillations.  PIC  is 
using  Ax  =  1 6Xp  with  10  cells  across  the  gap.  The  top  figure  in  a),  b)  and  c)  shows  the  initial 
distribution  (blue)  and  the  finial  distributions;  black-PIC,  red-BIT  and  magenta-BIT  corrected  PIC 
with  extended  cells  and  regulation.  The  parameters  are  a)  At  =  Ax/ vmax/ 1 0  and  d  =  Ax/ 20.  b) 
At  =  Ax/vmax/5  and  d  =  Ax/ 10  and  c)  At  =  Ax/vmax/2  and  d  —  Ar/10. 
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Figure  3:  Numerical  heating  test  with  4000  test  particles  uniformly  sampled  in  space  and  sampled 
from  a  Gaussian  in  velocity.  Black  line  -  PIC,  red  line  -  BIT,  magenta  line  -  BIT  Corrected  PIC.  a) 
At  =  Ax/vmax/10  and  d  =  Ax/20,  b)  At  =  Ax/vmax/5  and  d  —  Ax/ 10  and  c)  At  =  Ax/vmax/2  and 
d  =  Ax/10 


Collisional  Plasmas  (Christlieb,  Olson)  In  order  to  account  for  collisions,  the  right  hand  side 
of  equation  (1)  becomes 


5/ 

St 


L?  Is  vcdCldy*, 


where  v  is  the  relative  velocity  between  two  particles,  a  is  the  scattering  cross  section  and  d£t  is  the 
differential  solid  angle  through  which  particles  scatter.  The  two  product  terms  in  the  parentheses 
represent  a  source  and  loss  of  particles  which  scatter  into  and  out  of  a  given  velocity  through  colli¬ 
sional  events.  This  convolution  integral  over  velocity  space  is  expensive  to  compute.  One  common 
approach  to  reducing  the  cost  is  to  use  Monte  Carlo  techniques. 
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In  this  work  we  have  developed  an  innovative  grid-free  approach  to  Monte  Carlo  methods  for  gas 
dynamics  problems.  Again,  we  model  /  as  a  collection  of  macro  particles  in  phase  space.  The 
evolution  of  the  new  equation  (1)  is  carried  out  using  operator  splitting  in  time,  i.e.,  we  evolve  the 
macro  particles,  which  describe  the  transport  of  /  through  phase  space  and  then  correct  the  result¬ 
ing  distribution  functions  for  collisional  events  using  a  Monte  Carlo  method.  It  has  been  shown 
that  this  method  converges  to  the  solution  of  the  modified  equation  (1)  as  the  particle  number  goes 
to  infinity  .  Novel  aspects  in  this  work  are  that  particles  are  clustered  adaptively  using  octrees  at 
each  time  step  and  that  macro  flow  quantities  are  tracked  and  updated  in  the  nodes  of  the  tree  using 
grid-free  interpolation. 

The  method  has  been  applied  to  several  benchmark  problems  such  as  Couette  flow,  thermal  Cou- 
ette  flow,  and  flow  past  a  lifting  body.  In  addition,  we  have  applied  the  method  to  the  study  of  the 
formation  of  a  Bose-Einstein  condensate. 

Couette  and  Thermal  Couette  For  the  velocity  driven  Couette  flow,  we  compared  the  results  of 
our  grid-free  DSMC  with  a  thoroughly  tested  code,  MONACO.1  We  ran  several  test  cases  with 
different  Knudsen  numbers,  Kn.  The  two  approaches  are  in  excellent  agreement,  shown  in  figure 
4  (a)  and  (b).  Figure  4  (a)  is  the  temperature  profile  across  the  gap  for  the  velocity  driven  Couette 
flow  and  figure  4  (b)  is  the  velocity  profile  across  the  gap  for  the  velocity  driven  flow. 

Figure  4  (c)  shows  the  results  of  the  two  methods  for  the  temperature  driven  Couette  flow.  Again, 
the  results  of  the  two  methods  are  in  excellent  agreement.  The  results  for  flow  past  a  flat  plate  are 


Figure  4:  Simulation  results  for  grid  based  and  grid-free  Monte  Carlo  simulations  of  ther¬ 
mal  and  velocity  driven  Couette  flow.  The  images  show  results  for  a  range  of  densities,  kn  e 
{0.01,0.1, 1, 10, 100}.  Note  that  kn  is  the  mean  free  path  over  a  characteristic  length,  as  kn  —>  °° 
the  density  goes  to  zero.  On  the  far  right  the  plot  is  of  the  temperature  profile  for  the  thermal 
Couette  flow.  The  middle  and  left  plots  are  the  temperature  and  velocity  profiles  for  the  velocity 
driven  Couette  flow. 

also  in  good  agreement  with  previous  work.  Eventually  this  method,  in  combination  with  other 
grid-free  tools,  will  be  used  to  simulate  the  behavior  of  complex  multi-species  collisional  plasmas 
in  a  grid-free  framework.  Our  first  papers  on  this  method  is  complete  and  accepted  to  JCP  (May 
2008). 

'Grid  based  DSMC  Code  developed  by  Professor  Iain  D.  Boyd’s  research  group.  University  of  Michigan  Aerospace 
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Transitions 

AFRL/DEHE  is  collaborating  with  the  PI  on  transiting  BIT  corrected  PIC  into  AFRL-ICEPIC. 
New  Discoveries 

Numerical  Heating  in  Particle  codes  can  be  controlled  with  a  regularized  sub-cell  method  (PI). 
High  order  corrections  to  regularization  have  an  important  advantage  near  interfaces  (PI). 

High  order  corrections  have  an  important  advantage  near  interfaces  (PI). 

Analysis  of  oscillator  problem  points  to  why  regularization  is  useful  away  from  boundaries  (PI  and 
Dr.  Cartwright).  (Details  at  08  Review) 

High  order  correction  steps  inside  SDC  is  an  advantage  in  both  efficiency  and  accuracy.  (PI  and 
Dr.  Qiu)  (Details  at  08  Review) 

SDC  is  advantage  in  BIT  codes  even  when  using  point  insertion.  (PI  and  Dr.  Ong)  (Details  at  08 
Review) 
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